function qr_least_squares()
% QR_LEAST_SQUARES 基于QR分解的最小二乘求解
% 
% 演示使用QR分解求解最小二乘问题的方法和优势

fprintf('=== 基于QR分解的最小二乘求解 ===\n\n');

%% 1. QR分解求解最小二乘的原理
fprintf('1. QR分解求解最小二乘的原理\n');
fprintf('对于超定系统 Ax = b，其中 A = QR\n');
fprintf('最小二乘解: x = R^(-1) * Q^T * b\n\n');

% 示例问题
A = [1, 1;
     1, 2;
     1, 3;
     1, 4;
     1, 5];
b = [2.1; 2.9; 4.2; 4.8; 6.1];

fprintf('系数矩阵 A (%d×%d):\n', size(A));
disp(A);
fprintf('观测向量 b:\n');
disp(b');

% QR分解
[Q, R] = qr(A, 0);  % 经济型QR分解

fprintf('Q矩阵 (%d×%d):\n', size(Q));
disp(Q);
fprintf('R矩阵 (%d×%d):\n', size(R));
disp(R);

% 求解最小二乘
x_qr = R \ (Q' * b);
fprintf('QR分解解: x = [%.6f; %.6f]\n', x_qr);

% 验证解
residual = A * x_qr - b;
fprintf('残差范数: ||r|| = %.6f\n', norm(residual));

%% 2. 与正规方程法的比较
fprintf('\n2. 与正规方程法的比较\n');
fprintf('数值稳定性和精度对比\n\n');

% 正规方程解
AtA = A' * A;
Atb = A' * b;
x_normal = AtA \ Atb;

fprintf('正规方程解: x = [%.6f; %.6f]\n', x_normal);
fprintf('两种方法的差异: %e\n', norm(x_qr - x_normal));

% 条件数比较
fprintf('条件数比较:\n');
fprintf('cond(A) = %.2e\n', cond(A));
fprintf('cond(A^T*A) = %.2e\n', cond(AtA));
fprintf('cond(R) = %.2e\n', cond(R));

%% 3. 病态问题的处理
fprintf('\n3. 病态问题的处理\n');
fprintf('QR分解在病态问题上的优势\n\n');

% 构造病态矩阵
epsilon = 1e-10;
A_ill = [1, 1;
         1, 1+epsilon;
         1, 1+2*epsilon;
         1, 1+3*epsilon;
         1, 1+4*epsilon];
b_ill = [1; 2; 3; 4; 5];

fprintf('病态矩阵条件数: %.2e\n', cond(A_ill));

% QR分解求解
[Q_ill, R_ill] = qr(A_ill, 0);
x_qr_ill = R_ill \ (Q_ill' * b_ill);

% 正规方程求解
AtA_ill = A_ill' * A_ill;
Atb_ill = A_ill' * b_ill;
x_normal_ill = AtA_ill \ Atb_ill;

% MATLAB内置求解器
x_matlab = A_ill \ b_ill;

fprintf('病态问题求解结果:\n');
fprintf('QR分解解: x = [%.6f; %.6f]\n', x_qr_ill);
fprintf('正规方程解: x = [%.6f; %.6f]\n', x_normal_ill);
fprintf('MATLAB解: x = [%.6f; %.6f]\n', x_matlab);

% 残差比较
residual_qr = norm(A_ill * x_qr_ill - b_ill);
residual_normal = norm(A_ill * x_normal_ill - b_ill);
residual_matlab = norm(A_ill * x_matlab - b_ill);

fprintf('残差比较:\n');
fprintf('QR分解: %e\n', residual_qr);
fprintf('正规方程: %e\n', residual_normal);
fprintf('MATLAB: %e\n', residual_matlab);

%% 4. 秩亏问题
fprintf('\n4. 秩亏问题\n');
fprintf('当A不是列满秩时的处理\n\n');

% 构造秩亏矩阵
A_rank_def = [1, 2, 3;
              2, 4, 6;
              1, 2, 3.001;
              3, 6, 9];
b_rank_def = [1; 2; 1.5; 3];

fprintf('秩亏矩阵 A (%d×%d):\n', size(A_rank_def));
disp(A_rank_def);
fprintf('矩阵的秩: %d\n', rank(A_rank_def));
fprintf('理论秩: %d\n', size(A_rank_def, 2));

% 使用SVD处理秩亏问题
[U, S, V] = svd(A_rank_def, 0);
s = diag(S);
fprintf('奇异值: [%.2e, %.2e, %.2e]\n', s);

% 确定数值秩
tol = max(size(A_rank_def)) * eps(max(s));
rank_numerical = sum(s > tol);
fprintf('数值秩: %d (容差: %.2e)\n', rank_numerical, tol);

% 最小范数解
x_pinv = pinv(A_rank_def) * b_rank_def;
fprintf('最小范数解: x = [%.6f; %.6f; %.6f]\n', x_pinv);

%% 5. 增量QR分解
fprintf('\n5. 增量QR分解\n');
fprintf('当矩阵增加行或列时的高效更新\n\n');

% 基础矩阵
A_base = randn(6, 3);
[Q_base, R_base] = qr(A_base, 0);

fprintf('基础矩阵: %d×%d\n', size(A_base));

% 增加一行
new_row = randn(1, 3);
A_extended = [A_base; new_row];

fprintf('增加一行后: %d×%d\n', size(A_extended));

% 更新QR分解
[Q_updated, R_updated] = update_qr_add_row(Q_base, R_base, new_row);

% 验证更新
[Q_direct, R_direct] = qr(A_extended, 0);
fprintf('更新QR与直接QR的差异:\n');
fprintf('Q差异: %e\n', norm(Q_updated - Q_direct, 'fro'));
fprintf('R差异: %e\n', norm(R_updated - R_direct, 'fro'));

%% 6. 删除行的QR更新
fprintf('\n6. 删除行的QR更新\n');
fprintf('从QR分解中删除一行\n\n');

% 删除第3行
row_to_delete = 3;
A_reduced = A_extended;
A_reduced(row_to_delete, :) = [];

fprintf('删除第%d行后: %d×%d\n', row_to_delete, size(A_reduced));

% 更新QR分解
[Q_reduced, R_reduced] = update_qr_delete_row(Q_updated, R_updated, row_to_delete);

% 验证
[Q_direct_red, R_direct_red] = qr(A_reduced, 0);
fprintf('删除行更新的差异:\n');
fprintf('Q差异: %e\n', norm(Q_reduced - Q_direct_red, 'fro'));
fprintf('R差异: %e\n', norm(R_reduced - R_direct_red, 'fro'));

%% 7. 约束最小二乘
fprintf('\n7. 约束最小二乘\n');
fprintf('等式约束下的最小二乘问题\n\n');

% 问题: min ||Ax - b||² subject to Cx = d
A_const = randn(8, 4);
b_const = randn(8, 1);
C = randn(2, 4);  % 约束矩阵
d = randn(2, 1);  % 约束右端

fprintf('约束最小二乘问题:\n');
fprintf('A: %d×%d, C: %d×%d\n', size(A_const), size(C));

% 使用拉格朗日乘数法
% [A^T*A  C^T] [x]   [A^T*b]
% [C      0  ] [λ] = [d     ]

KKT_matrix = [A_const'*A_const, C'; C, zeros(size(C, 1))];
KKT_rhs = [A_const'*b_const; d];

solution = KKT_matrix \ KKT_rhs;
x_constrained = solution(1:size(A_const, 2));
lambda = solution(size(A_const, 2)+1:end);

fprintf('约束最小二乘解: x = [%.4f; %.4f; %.4f; %.4f]\n', x_constrained);
fprintf('拉格朗日乘数: λ = [%.4f; %.4f]\n', lambda);

% 验证约束满足
constraint_residual = C * x_constrained - d;
fprintf('约束残差: ||Cx - d|| = %e\n', norm(constraint_residual));

%% 8. 总最小二乘
fprintf('\n8. 总最小二乘\n');
fprintf('当系数矩阵A也有误差时的处理\n\n');

% 生成带噪声的数据
A_true = [1, 1; 1, 2; 1, 3; 1, 4];
x_true = [1; 2];
b_true = A_true * x_true;

% 添加噪声
noise_level = 0.1;
A_noisy = A_true + noise_level * randn(size(A_true));
b_noisy = b_true + noise_level * randn(size(b_true));

fprintf('噪声水平: %.2f\n', noise_level);
fprintf('真实参数: x = [%.1f; %.1f]\n', x_true);

% 普通最小二乘
x_ols = A_noisy \ b_noisy;
fprintf('普通最小二乘: x = [%.4f; %.4f]\n', x_ols);

% 总最小二乘 (使用SVD)
[U_tls, S_tls, V_tls] = svd([A_noisy, b_noisy]);
V_tls_22 = V_tls(end, end);
V_tls_12 = V_tls(1:end-1, end);
x_tls = -V_tls_12 / V_tls_22;

fprintf('总最小二乘: x = [%.4f; %.4f]\n', x_tls);

% 误差比较
error_ols = norm(x_ols - x_true);
error_tls = norm(x_tls - x_true);
fprintf('误差比较:\n');
fprintf('普通最小二乘误差: %.4f\n', error_ols);
fprintf('总最小二乘误差: %.4f\n', error_tls);

%% 9. 性能基准测试
fprintf('\n9. 性能基准测试\n');
fprintf('不同方法的计算时间比较\n\n');

sizes = [100, 200, 500, 1000];
fprintf('矩阵规模    QR分解    正规方程    SVD分解\n');
fprintf('--------    ------    --------    -------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 50);
    A_bench = randn(m, n);
    b_bench = randn(m, 1);
    
    % QR分解法
    tic;
    [Q_bench, R_bench] = qr(A_bench, 0);
    x_qr_bench = R_bench \ (Q_bench' * b_bench);
    time_qr = toc;
    
    % 正规方程法
    tic;
    AtA_bench = A_bench' * A_bench;
    Atb_bench = A_bench' * b_bench;
    x_normal_bench = AtA_bench \ Atb_bench;
    time_normal = toc;
    
    % SVD分解法
    tic;
    x_svd_bench = pinv(A_bench) * b_bench;
    time_svd = toc;
    
    fprintf('%6d      %6.4f    %8.4f    %7.4f\n', ...
            m, time_qr, time_normal, time_svd);
end

%% 10. 应用实例：信号去噪
fprintf('\n10. 应用实例：信号去噪\n');
fprintf('使用QR分解进行信号的多项式去噪\n\n');

% 生成含噪信号
t = linspace(0, 1, 100)';
signal_true = sin(2*pi*t) + 0.5*cos(6*pi*t);
noise = 0.3 * randn(size(t));
signal_noisy = signal_true + noise;

% 构造多项式基函数矩阵
poly_degree = 10;
A_poly = zeros(length(t), poly_degree + 1);
for j = 0:poly_degree
    A_poly(:, j+1) = t.^j;
end

fprintf('信号长度: %d\n', length(t));
fprintf('多项式次数: %d\n', poly_degree);
fprintf('信噪比: %.2f dB\n', 20*log10(std(signal_true)/std(noise)));

% QR分解去噪
[Q_denoise, R_denoise] = qr(A_poly, 0);
coeffs_denoise = R_denoise \ (Q_denoise' * signal_noisy);
signal_denoised = A_poly * coeffs_denoise;

% 计算去噪效果
mse_original = mean((signal_noisy - signal_true).^2);
mse_denoised = mean((signal_denoised - signal_true).^2);
improvement = 10 * log10(mse_original / mse_denoised);

fprintf('去噪效果:\n');
fprintf('原始MSE: %.6f\n', mse_original);
fprintf('去噪后MSE: %.6f\n', mse_denoised);
fprintf('改善: %.2f dB\n', improvement);

% 绘制结果
figure('Name', '信号去噪 - QR分解');
plot(t, signal_true, 'g-', 'LineWidth', 2);
hold on;
plot(t, signal_noisy, 'r.', 'MarkerSize', 8);
plot(t, signal_denoised, 'b-', 'LineWidth', 2);
xlabel('时间');
ylabel('幅度');
title('基于QR分解的信号去噪');
legend('真实信号', '含噪信号', '去噪信号', 'Location', 'best');
grid on;

end

function [Q_new, R_new] = update_qr_add_row(Q, R, new_row)
% UPDATE_QR_ADD_ROW 向QR分解中添加一行
% 
% 输入: Q, R - 原始QR分解, new_row - 新增行
% 输出: Q_new, R_new - 更新后的QR分解

[m, n] = size(R);

% 扩展Q矩阵
Q_extended = [Q; zeros(1, n)];

% 计算新行在当前列空间中的投影
projection = new_row * Q;
residual = new_row - projection;
residual_norm = norm(residual);

if residual_norm > eps
    % 添加新的正交向量
    q_new = residual / residual_norm;
    Q_new = [Q_extended, [zeros(m, 1); 1]];
    
    % 更新R矩阵
    R_new = [R; projection];
else
    % 新行在列空间中，不增加维数
    Q_new = Q_extended;
    R_new = [R; projection];
end

end

function [Q_new, R_new] = update_qr_delete_row(Q, R, row_index)
% UPDATE_QR_DELETE_ROW 从QR分解中删除一行
% 
% 输入: Q, R - 原始QR分解, row_index - 要删除的行索引
% 输出: Q_new, R_new - 更新后的QR分解

% 删除Q中的对应行
Q_reduced = Q;
Q_reduced(row_index, :) = [];

% 使用Givens旋转恢复上三角形式
[m_new, n] = size(Q_reduced);
R_new = R;

% 重新正交化（简化版本）
[Q_new, R_new] = qr(Q_reduced * R, 0);

end